load("data/ocdist_16.Rdata")
# plot(pg)
dat <- pg@data

#abramson data
library(haven)
library(maptools)
library(sf)
library(rgdal)
library(sp)

#Euroatlas shapes
#eur12 <- readOGR("Shapefiles/Euroatlas/1200/sovereign_states.shp", layer="sovereign_states")
#proj4string(eur12) <-  "+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"
#proj4string(eur16)

lee <- read.csv2("data/Parliaments_lee_v4.csv")

tab <- read_dta("data/table_6.dta")

tab <- subset(tab, select=c(pid, Year, lnparl))

#merging LEE with Abramson
colnames(tab)[2] <- "year"
lee <- cbind(lee, tab)


#lat
lee$lat <- as.character(lee$lat)

lee$lat <- as.numeric(gsub(",", ".", lee$lat))
lee$lat <- as.numeric(as.character(lee$lat))


#long
lee$long <- as.character(lee$long)

lee$long <- as.numeric(gsub(",", ".", lee$long))
lee$long <- as.numeric(as.character(lee$long))

lee2 <- aggregate(lee[4:5], by=list(lee$Name), FUN=max, na.rm=T)
colnames(lee2)[1] <- "Name"

lee <- subset(lee, select=c(pid, Name, year, lnparl))
lee <- merge(lee, lee2, by=c("Name"), all.x=T)

#####SETTING THE PROJECTION
lee2 <- lee


lee <- lee[is.infinite(lee$lat)==F,]
lee <- lee[is.infinite(lee$long)==F,]

coordinates(lee) <- cbind(lee$long, lee$lat)

## tiff('tmp.tiff')
## plot(lee, add=T, col="red")
## dev.off()

#now I need to get the gid for each year

proj4string(eur16)

proj4string(lee) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"

#taking the parliaments for a grid cell
l <- over(lee, pg)

lee <- cbind(lee@data, l$gid)
colnames(lee)[7] <- "gid"

#Euroatlas polygons over grid cells
l <- over(pg, eur16)

l <- subset(l, select=c(owner_id))

###WORLD SHAPES
library(cshapes)

cs <- cshp(date=as.Date("1950-6-30"), useGW=T)

proj4string(cs)
proj4string(pg)
cs <- spTransform(cs, "+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0")
cs <- subset(cs, select=c(GWCODE))

l2 <- over(pg, cs)
l2$GWCODE[is.na(l2$GWCODE)==T] <- 999

#merging on population and river distance
pop <- read_dta("data/hydepop_1900_1940.dta")

pop <- subset(pop, year==1900)

pg <- pg@data

pg <- merge(pg, pop, by=c("gid"))

#mergin on euroatlas
pg <- cbind(pg, l)

#merging on cshapes
pg <-cbind(pg, l2)
#merging on ocean distance vars
suit <- read_dta("dist2suitable.dta")

pg <- merge(pg, suit, by=c("gid"))


#merging on new predictive vars
pred <- read_dta("predportout.dta")

pg <- merge(pg, pred, by=c("gid"))



#splitting up for each year and merging on parliaments
pg1 <- pg
pg1$year <- 1200

pg2 <- pg
pg2$year <- 1300

pg3 <- pg
pg3$year <- 1400

pg4 <- pg
pg4$year <- 1500

pg5 <- pg
pg5$year <- 1600

pg6 <- pg
pg6$year <- 1700


pg <- rbind(pg1, pg2, pg3, pg4, pg5, pg6)

lee$parl <- 1
head(lee)

lee <- subset(lee, select=c(gid, year, parl,lnparl ))

pg <- merge(pg, lee, by=c("gid", "year"), all.x=T)

pg$parl[is.na(pg$parl)==T] <- 0
pg$lnparl[is.na(pg$lnparl)==T] <- 0

###Operationalizing
table(pg$parl)

#distance
pg$two_port_all_pred1 <- log(pg$two_port_all_pred1/1000+0.00001)

pg$two_port_all_pred1  <- pg$two_port_all_pred1/100


pg$dist2suitable <- log(pg$dist2suitable/1000+0.00001)

#population
pg$lpop <- log(pg$sum+0.0001)


head(pg)


#port suitable (parliament binary)
summary(mod1 <- lm(parl ~ two_port_all_pred1 + as.factor(year)-1 , data=pg))
summary(mod1,cluster = c("gid"))


summary(mod2 <- lm(parl ~ two_port_all_pred1 + lpop + city + as.factor(year)-1, data=pg))
summary(mod2,cluster = c("gid"))


#save clustered se's
cluster_se1 <- as.vector(summary(mod1,cluster = c("class_id"))$coefficients[,"Std. Error"])
cluster_se2 <- as.vector(summary(mod2,cluster = c("class_id"))$coefficients[,"Std. Error"])


library(stargazer)
library(sjPlot)
# print stargazer output with robust standard errors
stargazer(mod1,mod2,
          type = "text",se = list(cluster_se1, cluster_se2,
                                  cluster_se3), out="models1.xls")
?stargazer

####YEARS
summary(mod1 <- lm(parl ~ two_port_all_pred1 + lpop + city , data=pg[pg$year==1200,]))
summary(mod2 <- lm(parl ~ two_port_all_pred1 + lpop + city , data=pg[pg$year==1300,]))
summary(mod3 <- lm(parl ~ two_port_all_pred1 + lpop + city , data=pg[pg$year==1400,]))
summary(mod4 <- lm(parl ~ two_port_all_pred1 + lpop + city , data=pg[pg$year==1500,]))
summary(mod5 <- lm(parl ~ two_port_all_pred1 + lpop + city , data=pg[pg$year==1600,]))
summary(mod6 <- lm(parl ~ two_port_all_pred1 + lpop + city , data=pg[pg$year==1700,]))




